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Abstract 

Phase-Rectified Signal Averaging (PRSA) was shown to be a powerful tool for the 
study of quasi-periodic oscillations and nonlinear effects in non-stationary signals. 
Here we present a bivariate PRSA technique for the study of the inter-relationship 
between two simultaneous data recordings. Its performance is compared with tra- 
ditional cross-correlation analysis, which, however, does not work well for non- 
stationary data and cannot distinguish the coupling directions in complex nonlinear 
situations. We show that bivariate PRSA allows the analysis of events in one signal 
at times where the other signal is in a certain phase or state; it is stable in the 
presence of noise and impassible to non-stationarities. 
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1 Introduction 

Many natural systems generate periodicities on different time scales because 
some of their components form closed regulation loops in addition to causal 
linear control chains. In biology and physiology, cardio-respiratory rhythms, 
rhythmic motions of limbs in walking, rhythms underlying the release of hor- 
mones and gene expression, membrane potential oscillations, oscillations in 
neuronal signals, and circadian rhythms are just a few examples (see, e.g., 
p1l2] ). Oscillations also occur in geophysical data, e.g., for the El- Nino phe- 
nomenon, sunspot numbers, and ice age periods [3j. In many cases several 
signals from different components of the complex system can be recorded si- 
multaneously. For understanding the control chains and loops in the system. 
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we want to know how periodicities in the signals are generated by (possibly di- 
rected and/or nonlinear) interactions between its components. Consequently, 
there is a need for identifying periodicities in one recorded signal together with 
the direction of causal relations to periodicities in other signals. 

Cross-correlation analysis and transfer function analysis are traditional tools 
for this type of analysis. However, there are three major drawbacks of these 
methods: (i) only rather stationary data can be studied, (ii) a linear rela- 
tionship between the signals is usually assumed, and (iii) the identification of 
causalities is hindered by the fact that the exchange of the two signals under 
study is identical with time inversion. We thus propose a method which helps 
to overcome these problems. 

Non-stationarities are a major problem in the analysis of signals recorded from 
complex systems over a prolonged period of time [H[5|5f7f5] . Many internal and 
external perturbations are continuously influencing the system and causing 
interruptions of the periodic behavior. The interruptions often 'reset' the reg- 
ulatory mechanisms resulting in phase de-synchronization of the oscillations. 
The signals thus become quasi-periodic, consisting of many periodic patches as 
well as noise and trends. Cross-correlation and transfer function techniques are 
thus problematic. In addition, there might be causal inter-relations between 
two signals that cannot be revealed by these methods. For illustration, let us 
assume that a large increase and a large decrease in signal X (trigger signal) 
cause the same specific effect in signal Y (target signal), while there is no 
such effect in F if X remains unchanged. In this situation with an essentially 
nonlinear coupling between the signals, both, cross-correlation analysis and 
spectral analysis cannot reveal the effect. They show the superposition of the 
two branches of the interaction with opposite signs, i.e., no effect. Even if the 
effects on signal Y were different for increases and decreases of signal X, one 
could see some relation but could not distinguish the two effects. Hence, one 
needs a method that can separately study effects in signal Y which might oc- 
cur in response to different causes in signal X, and vice versa. A separation of 
effects with different typical duration or frequency scale seems also appropri- 
ate for distinguishing frequency-band selective inter-relations between signals 
X and Y. 

Our approach for extracting inter-relations between two or more simultaneous 
data recordings from a complex system is based on the phase-rectified signal 
averaging technique (PRSA) piTO] . which was shown to be a powerful tool 
for the study of quasi-periodic oscillations in noisy, non-stationary signals. 
The original method extracts the features in one signal before and after in- 
creases in the same signal (or, alternatively, decreases). This way, information 
on characteristic quasi-periodicities, short-term correlations, and time inver- 
sion asymmetry (causality) is extracted, while non-stationarities and noise are 
eliminated. The advanced approach introduced here extracts the features in 
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one signal before and after increases in another signal. Thus, the inter-relation 
between both signals can be studied separately for both coupling directions, 
both time directions and independent of non-stationarities and noise. 

The paper is organized as follows. In Section 2 we describe both, the univariate 
and the bivariate PRSA method. Section 3 is dedicated to the comparison of 
the bivariate PRSA method with the traditional cross-correlation analysis. We 
also address pitfalls and drawbacks of the cross-correlation analysis that are 
often overlooked. In Section 4 we discuss three model examples and quantify 
the capacity of the bivariate PRSA method for the detection of nonlinear in- 
teractions and quasi-periodicities. Finally, we summarize and discuss possible 
applications in Section 5. 



2 PRSA methods 

2.1 Univariate PRSA 

Let X = [xi), i = 1, . . . , N he a. long time series representing the signal under 
investigation. In addition to periodicities and correlations of interest, X may 
contain non-stationarities, noise and recording artifacts. One example for such 
a signal is the series of time intervals between successive heartbeats determined 
from a long-term EGG (electrocardiogram) of a patient in a hospital. Since 
the most pronounced peak in the EGG used for heartbeat interval determina- 
tion is called the R-peak, these time series are often denoted as RR-interval 
(RRI) time series. Phase Rectified Signal Averaging was shown to reduce the 
signal to a much shorter sequence keeping all relevant quasi-periodicities but 
eliminating non-stationarities, artifacts, and noise [9]. The PRSA algorithm 
consists of three major steps as illustrated in Fig. 1. 

Step 1. Anchor points in the time series are defined according to specific fea- 
tures, e.g., increases (or, alternatively, decreases) in the time series (see Fig. 
1(a)). I.e., a point Xj qualifies as an anchor point if 



for triggering on increases or decreases, respectively. In order to study a lower 
frequency regime, averages of T successive values of Xi are compared Typi- 
cally half of all points of the time series will qualify as anchor points. In general, 
quasi-periodic oscillations in a noisy time series X will result in anchor points 
predominantly found in the phase of the steepest ascent (or decent for the 
second alternative in Eq. ([1])), i.e., when the phase of the signal itself is close 
to (or close to tt). The phase information of the oscillations is thus obtained 
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Fig. 1. Illustration of PRSA technique: (a) Anchor points selected in the time series 
(xj), here: increases; only the first 5 anchors are marked, (b) Windows (surroundings, 
here: only first 4 shown) of length 2L (here: L = 16) defined around each anchor 
point, (c) Surroundings of all anchor points on top of each other (here: only first 25 
shown), (d) PRSA curve x{k) from averaging over all surroundings; the parameter 
L is increased to L = 32 in order to improve the visibility of the slow periodicity. 
The original signal is 1// noise (generated with the Fourier filtering method) with 
two additional quasi-periodicities with characteristic frequencies / = 0.05/ At and 
/ = 0.3/ At; phase jumps are inserted after an average number of four periods (from 
i). 



from the signal itself, and the signal can be phase-rectified using the anchor 
points. Note, that in principle any boolean valued function may be used for 
the definition of anchor points, where true is associated with an anchor while 
false is not. This allows studying more complex structures in signals. 

Step 2. Windows, i.e., surroundings, of length 2L around each anchor point 
Xi^, V = 1, . . . , M, are identified (see Fig. [11(b)); M is the total number of 
regarded anchor points. The surrounding of is 

^i„-L^ ^i^-L+ly • • • ; • • • 5 2;j^+i_2, (2) 

The parameter L has to be chosen larger than the expected coherence time 
of the periodicities in the signal; it must definitely exceed the period of the 
slowest oscillation that one wants to detect. All anchor points with indices v 
smaller than L + 1 and larger than A^ — L + i.e., at the very beginning and at 
the end of the time series, have incomplete surroundings. The same holds for 
windows containing missing data points due to, e.g., measurement artifacts, 
instrument failure, or outliers. 
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Step 3. All windows z/, v = 1, . . . , M are aligned at their anchor points Xi^, 
and the phase-rectified signal average x{k) is obtained by averaging over all 
windows (see Figs. 1(c) and (d)), 

1 M 

PRSAx(A;) = x(A;) = — k = -L, . . . ,0, . . . , L - 1. (3) 

If Xi^^k is a missing data point, it is replaced by 0, and M is substituted by 
denoting the number of non-missing points at position k. Including windows 
with missing data points yields better statistics and allows investigation of 
time series with a few artifacts. In general a well-behaved average x{k) can be 
expected when there are at least 100 to 1000 anchor points, i.e., = 200 to 
A^ = 2000 for the length of the record. 

In the average ([3]), non-periodic components (not phase-synchronized with 
the anchor points), i.e., non-stationarities, non-identified artifacts, and noise, 
cancel out. Only events that have a fixed phase relationship with the anchor 
points, i.e., periodicities and quasi-periodicities, 'survive' the procedure (see 
Fig. [T]). The PRSA signal x{k) represents the most important features of the 
original data containing all quasi-periodicities aligned with phase zero in the 
center (at k = 0). Applying the PRSA before traditional spectral analysis 
significantly improves the quality of the spectra in the presence of noise and 
non-stationarities piTT] . Differences between PRSA curves obtained by ap- 
plying either of the two criteria in Eq. ([1]) will indicate missing time reversal 
symmetry of the original signal. Hence, nonlinear and non time-reversal in- 
variant processes, with different phenomena occurring during increasing and 
decreasing parts, can be studied in detail. Optionally, it might be meaning- 
ful to weight the windows according to some criteria, e.g., according to the 
magnitude of changes at anchor positions in the trigger signal. With anchors 
defined at increases Eq. becomes 

M 

FRSAxik)=xik) = Y.Ci^x^^+k, k = -L,...,0,...,L-l. (4) 

1^=1 

with weights q^, e.g., = (xj^ - Xi^-i) / J2ff=i{xi^ - Xj^-i). Of course, other 
weights could be defined as well. 

2.2 Bivanate PRSA (BPRSA) 

Now, we suggest a generalization of the univariate PRSA for studying the 
inter-relations between two signals X and Y . If many signals are recorded 
simultaneously, representing the dynamics of the complex system, each pair 
can be characterized accordingly. 
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Fig. 2. Illustration of the BPRSA technique: According to Eq. ([T]) anchor points are 
(a) selected in a trigger signal X and (b) transferred to the target signal Y. Here only 
parts of much longer blood pressure (a) and heartbeat interval (b) recordings are 
shown. After averaging the windows around each anchor point in Y (red) according 
to Eq. ([5]), BPRSAx^y is obtained (c). Likewise, changes in X caused by increases 
in Y can be studied by exchanging trigger and target signal, BPRSAy^x in (d). 



The method is nearly identical with the univariate approach described in the 
previous subsection, except for the usage of different signals in step one (xj = 
trigger signal X) and in steps two and three (xj = target signal Y). Specifi- 
cally, anchor points i^, v = 1, . . . , M are defined for increases (or alternatively 
decreases) in the trigger signal X, i.e. (xj) (step 1), while surroundings are 
defined (step 2) and averaged (step 3) for the target signal y, i.e. {yi). This 
yields the bivariate phase rectified signal average y{k): 

T M 

BPRSAx-.y (A;) = y{k) = tj E = -L, . . . , 0, . . . , L - 1. (5) 

The transfer of the anchor points is illustrated in Figs. [2](a),(b). 

BPRSA is a non-symmetric algorithm, i.e., the exchange of trigger signal X 
and target signal Y will result in a different BPRSA curve, see |2](c),(d). More 
complex boolean or weighted anchor criteria, even ones based on more than one 
trigger signal, are possible. For example, the typical behavior of a target signal 
Y can be studied around points of time (anchors) with increases of signal Xi 
and positive values of signal X2. For a specific example of a conditional anchor 
criterion consider the three signals heartbeat intervals, respiratory phase, and 
blood pressure. One can study characteristic heartbeat intervals (target signal 
Y) around increasing systolic blood pressure (first trigger signal Xi) and at a 
certain respiratory phase (second trigger signal X2). First attempts revealed 
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promising results for the investigation of baroreflex properties and will be 
reported in a medical journal. 



3 BPRSA and cross-correlation analysis 



3.1 Cross- correlation analysis 



Although cross-correlation analysis is considered as a well established tool for 
the study of inter-relations between two signals in many applications, only a 
few authors have specifically studied its reliability [T2|13fl4| . For two discretely 
measured signals (xj) and i = 1, . . . ,N, the normalized cross-correlation 
function is most commonly defined as 

N-k 

pxy{k) = — V (xi - fXx){yi+k - fJ'y) for /c = 0, 1, . . . and (6a) 

Na^ay .^^ 

1 ^ 

Pxy{k) = — V (xi - - /i?;) for = -1, -2, . . . . (6b) 

Na^ay 



1/2 



Here, Pa = jj Hi=i and (Ta = J^fLiic^i ~ PaY are mean and standard 
deviation of both series a = x,y, respectively. This definition assumes that 
both and cTq, do not vary in time, i.e., they do not depend on the segments 
of the time series selected for the study. This corresponds to the assumption 
of weak stationarity. Strong stationarity additionally requires constancy of all 
other moments. For studies discussing the replacement of fix and fiy by local 
estimates, e.g. running averages, see [T51IT6] . Note, however, that some cross- 
correlations might be reduced or eliminated by this so-called pre-whitening 
procedure, which is therefore unsafe. 

Another problem of cross-correlation functions is that the exchange of the two 
signals X and Y corresponds to replacing k by —k, i.e., time inversion. Hence, 
causality relations between the two series can hardly be assessed. In general, 
the points of Pxy{k) are highly auto-correlated, e.g., pxy{k) is strongly corre- 
lated with Pxy{k + 1). I.e., neighboring points in pxy{k) are stronger correlated 
with each other than neighboring points in the original time series P^SfTT] . This 
self-correlation causes long living trends in pxy{k), e.g., a slow decay after a 
peak, which is at risk of misinterpretation. 

Furthermore, the sum in Eqs. runs over N — k terms, while it is divided 
by N instead of N — k. This procedure corresponds to a standard averag- 
ing procedure only in the limit of very long data {N oo). Nevertheless, 
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most statistical toolkits employ the definition ([6]), because the convolution 
theorem and fast Fourier transform can be used to speed up the calculations 
significantly in this case by application of the Wiener- Khinchin theorem. Some 
authors even argue for an increase in precision because the normalization 
reduces the mean-square variance of Pxy{k) (see, e.g., p2]). However, this non- 
matching prefactor results in a bias towards zero with increasing time lag k 
for small N, causing a triangular-shaped behavior of pxy{k). Consequently, the 
value of I A; I > for the center of a peak in Pxy{k) is systematically underes- 
timated |l13jJ. We are convinced that the correction factor N/{N — k) which 
transforms Pxy{k) from Eqs. (jnD into the correctly normalized cross-correlation 
function 

2 N-k 

CCFx,Y{k) = TTT 7T J2 i^i " /^x)(l/i+fc " Py) , (7) 

{N-k)axay fr[ 

must be used to obtain reliable results except for very long data. 

If the considered data is not fully stationary, one might want to use only the 
values Xi with i = 1, . . . ,N — k and yi with z = + 1, . . . , for calculating px, 
Py, ax, and ay. This approach is known as local cross- correlation in literature; 
it is equivalent to the Pearson rxy (product-moment) correlation coefficient 
for the two overlapping pieces. Since the partial means and standard devia- 
tions will depend on k, the computational effort is significantly increased. The 
bias mentioned in the previous paragraph is not completely removed in this 
approach [13] (although it is weaker than for the standard definitions ([6])). 
In addition, problems with autoregressive moving average processes (ARMA) 
were reported [17j. Since the cross-correlation approach does not work well for 
non- stationary data anyway, we do not consider local cross-correlation here. 



3.2 Interpretation of BPRSA curves 



In BPRSA, anchor points usually occur in all parts of the trigger signal X. 
The average of BPRSAx^y (fc) = y{k) for all k will thus be approximately the 
global average of the whole signal, i.e., Py. Consequently, subtraction of this 
mean from y{k) yields positive and negative values as in the cross- correlation 
function. y{k) thus be interpreted in a similar way as an unnormalized 

cross-correlation function. If one divides by the global standard deviation, cr^, 
the resulting quantity 

BPRSA5!:^)(A;) = ^^^^ ~ (8) 

ay 

is also normalized. It can be compared with CCF xyik) in Eq. ([7]) and inter- 
preted in a similar way. Note that - different from cross-correlation analysis 
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- this rescaling is just the last step, and idy does not enter into the calcu- 
lation of BPRSAx^y(fc)- Hence, the shape of the curve cannot be affected 
by non-stationarities, i.e., inaccurate /ly. There is no practical advantage of 
normalized BPRSA, unless the behavior of the curves for different signals, 
e.g., triggering directions {X — > Y) and (y — > X), shall be directly compared. 
However, the global mean /jy and global standard deviation ay might not exist 
due to non-stationarities and in this case normalization is not recommended. 

In some apphcations it is even preferred to study the unnormalized BPRSA 
curves. For example, in quantifying the action of blood pressure upon heart- 
beat regulation via the baroreflex mechanism in the human cardiovascular 
system, the variation of the time intervals between successive heartbeats in 
reaction to increases in blood pressure needs to be measured. In this case 
the units of both signals have to be kept, and the measure characterizing the 
baroreflex must have the unit ms/mmHg, i.e., time difference divided by pres- 
sure difference. In fact, cross-correlation studies can only yield either quantities 
without units (if normalized) or quantities which are products of both orig- 
inal units. Quantities with the unit of only one original series or their ratio 
(as needed for the baroreflex) cannot be obtained. Hence, there is no way to 
obtain a meaningful measure for the baroreflex from a cross-correlation analy- 
sis, although the baroreflex is a typical example of a meaningful inter-relation 
between two components of a complex system. 

Effects occurring in BPRSAx^y(/E) for /c > can be easily recognized as 
consequences of the triggering events in the trigger signal X. On the other 
hand, effects seen in BPRSAx^y(/c) for k < are likely to be causes for 
the actual triggering events. Note that a similar conclusion is also valid for 
the cross-correlation function CCF x,Y{k), since effects observed for k > 
and k < are probably due to interactions from signal X onto signal Y and 
vice versa. However, BPRSA allows separating these causality effects from 
nonlinear effects, as we will see in the following. 

Altogether, four BPRSA curves can be defined, compared with one 
cross-correlation function: BPRSAj^_^y(/c) (triggering on increases in X), 
BPRSA^_^y(fc) (triggering on decreases in X), BPRSA^_^^(A;), and 
BPRSAy>_^^(A;) (triggering on Y). By comparing these curves, additional in- 
formation on the linearity of the interactions and time reversal symmetry can 
be obtained. In the following we will use the symbols and \ for BPRSA 
curves obtained by triggering on increases and decreases in the trigger sig- 
nal only when necessary for distinction. In all other cases BPRSAx->y means 

If the interaction from signal X to signal Y is linear, we will find 
BPRSAj^_^y(A;) = — BPRSA^^y(/c), since increases and decreases in X must 
cause opposite effects in Y. Accordingly, BPRSA^_^j(-(/c) — — BPRSAyL^^^ (A;) 
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shows that the interaction from y to X is hnear. If the interaction between 
both signals is fully symmetric, time inversion is equivalent with exchanging 
the signals, BPRSA5^™(A;) = BPRSA^i:°™)(-^) and BPRSA>1:°™)(A;) = 
BPRSAyi^^™^(— fc). Deviations from this behavior show non-symmetric cou- 
pling as do deviations from CCFxxi^) = CCFx,y(— fc) in cross-correlation 
analysis. However, this can be checked independent of the linear or nonlin- 
ear character of the interactions between the signals. Note that normalized 
BPRSA must be considered in this case, Eq. ([8]). It is straightforward to 
write down similar relations for testing further hypothesis regarding the inter- 
relations between both signals. 

3.3 Comparison of cross- correlation analysis and BPRSA 

In this subsection we will see how the BPRSA overcomes the disadvantages 
of cross-correlation analysis described before. 

1. Causality and nonlinear interactions. As we have shown in the previous 
subsection, more information on the linearity or nonlinearity of the interac- 
tions and on time-reversal symmetry can be obtained from BPRSA curves 
than from the cross-correlation function. 

2. Time delays. The estimation of (positively or negatively) time-delayed inter- 
relations between both signals is straightforward, just as in cross-correlation 
analysis. 

3. Missing data and outliers. BPRSA can easily cope with missing data (e.g., 
measurement artifacts, instrument failure, or outliers) in both series X and 
Y. Invalid values in X just cannot become anchor points. Invalid values in Y 
will be disregarded (see text following Eq. ([3])). 

4. (Non-)stationarity of the data. In the definition of BPRSA (Section 2, in 
particular Eq. ([5])) neither means nor standard deviations of both signals X 
and Y are needed. Hence, no direct problems arise for non-stationary data. 
In particular data with a piecewise constant trend, which is often observed in 
medical data recordings, will cause no problems in BPRSA, because Eq. ^ is 
a simple linear arithmetic averaging procedure. The deviations from a small 
or large local average will have the same weight in this averaging procedure. 
Hence, BPRSA does not need pre-whitening of the data before analysis. Cross- 
correlation analysis, on the other hand, will be disturbed severely by a piece- 
wise constant trend, because the deviations Xi — fi^ from the global average 
will be dominated by this trend (see Subsection 4.3 for an example). The same 
holds for an oscillating trend in the target signal Y which is uncorrelated with 
the trigger signal X. However, such a trend in X will selectively cause anchor 
points and thus disturb also BPRSA; consequently more anchor points, i.e.. 
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longer data, will be needed! 



A slowly varying, monotonous (e.g., polynomial) trend in the target signal will 
bend the BPRSA curve, since the local means are different in the beginning 
and at the end of the signal and in the beginning and at the end of each seg- 
ment. However, this bending is definitely not stronger than a similar bending 
of the cross-correlation function. Trends in the trigger signal X will modify 
the fraction of anchor points for increases and decreases, which has little effect 
on BPRSAx^y(A;) unless these trends are very strong. 

5. Enhanced auto- correlations. Unlike the cross-correlation function [T3p7] . 
which is often dominated by low frequencies, BPRSA does not show artificially 
enhanced auto-correlations. On the contrary, low frequencies are reduced due 
to the filtering characteristics (see next point). This makes BPRSA particu- 
larly attractive for studying signals with underlying 1//- rather than white 
noise. Note that l//-noise is prevalent, e.g., in medical and geophysical data. 

6. Filtering characteristics. Figure [3] compares the spectral properties of both, 
cross-correlation analysis and BPRSA. Since many interesting data contain 
long-term auto-correlations and are characterized by 1/ /-noise in their power 
spectra, P{f) ~ with /3 around 1, we start with two such noise series 
(see Fig. [a(a,b)) with ^ 0.7 and py ^ 1.3 (see Fig. (Sj^d)). The power 
spectrum of the cross-correlation function decays as (see Fig. [3](e)). It 
is thus dominated by low-frequency components. The BPRSA curve, on the 
other hand, yields a nearly flat power spectrum (see also Fig.[3](e)). Therefore, 
additional peaks and quasi-periodicities can be noticed and determined much 
easier. 

The filtering characteristics of BPRSA can be motivated as follows. The scaling 
behavior of the BPRSA spectrum is influenced by the anchoring procedure 
in the trigger signal and by the averaging of the target signal. We want to 
estimate the probability p{f) that an oscillating component with frequency 
/, Uf = Ay sm{27i ft) in the target signal Y affects BPRSAx-^y (A;) under the 
condition that an oscillation with the same frequency /, Xf = AxSm{2'n'ft) 
causes anchor points in the trigger signal X. Firstly, Xf has to cause anchor 
points at positions t^, meaning Xf{t^) has to be larger than Xj(t^ — At) ^ 
Xfity) — Atx'jr = Xfity) — At27ifAxCos{2nft„) for anchor criterion Eq. ([1^). 
This is a valid approximation except for very high frequencies /. The deviation 
Xf{t^) — Xfity — At) = A.t2'K f A^ cos{2ti ft becomes maximal for t^ = n/ f 
with any integer n. Since anchor points t^ are primarily generated at or close 
to phase zero of the considered component Xf, the later averaging is phase- 
rectifying in terms of the trigger signal. The value of the maxima Xf{ty) — 
Xf{tv — At) is 271 At f Ax and thus the probability px to anchor is proportional 
to Axf. On the other hand, the component yf has an effect proportional to its 
amplitude Ay due to the averaging procedure of Eq. and therefore Py Ay. 
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Fig. 3. Filter-properties of BPRSA and cross-correlation function for differently 
correlated noises with Hx = fJ-y = 0, ax = cry = 1. The time domain (a,b) and 
frequency domain (d) plots of l//^-noise with spectral exponents Px = 0.7 (black 
curves, shifted) and /3y = 1.3 (magenta curves) are shown next to the BPRSAx^y 
(black) and CCFx,y (blue) in (c) and their correctly normalized and correspond- 
ingly color-coded power spectra with /JsPRSAx-.y ~ (shifted) and PccFx y ~ 2 
in (e). The two long-term correlated noises were generated by Fourier filtering with 
different /?, starting both procedures with the same original white noise (not shown). 
For all spectra logarithmic binning and linear fitting (yellow dots and lines) were 
applied to estimate /3. 

The amplitude of the considered spectral components in BPRSAx-+y(/i;) is 
thus determined by A^Ayf. If we consider two signals X and Y consisting of 
correlated noise with power spectra 

Px{f)^Al^r^- and Pyif) ^ Al f-^- (9) 

we obtain 

PbPRSa(/) ~ {p.Pyf ~ AlfAl ~ ^ ^-;3BPaSA (^Q) 



with /?BPRSA = /3x + /3y — 2, yielding /3bprsa ~ if both (3x and /3y are close 
to one or their average is close to one. 
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4 Three illustrative examples 



Since BPRSA has significant advantages over cross-correlation analysis for 
studying data with 1// noise and/or nonlinear interaction as well as non- 
stationary data, one can imagine several applications. Here, we describe three 
specific situations and illustrate the performance of BPRSA on model data. 

4-1 White noises with linear relation 

We consider two independent white noise signals X = (xj) and Y = {tji) with 
zero mean and unit variance. Based on Y we generate the signal Y = {i/i) 
by introducing a linear unidirectional coupling with X in a certain frequency 
band. This is generated by calculating the linear combination of Y and one or 
more bandpass filtered components of X, 

y^ = y^ + J2^^M'\x). (11) 
j 

The bandpass filtering is done in Fourier space, and BPp\x) denotes the 
i-th element of the series obtained from the related j-th bandpass filter op- 
erator acting on X. The prefactors cj include the coupling strengths \cj\ and 
directions sgn(cj). Finally, Y is normalized to obtain zero mean and unit vari- 
ance. Fig. m^a) illustrates the original noise X, while Figs. Il](b,c) show Y and 
Z = (zi) for two different values of Ci and Cj = 0, Vj > 1. 

Different coupling strengths |ci| are reflected by different amplitudes of 
BPRSAx^a(fc) and a = Y, Z, while a different coupling direction results in a 
different sign of BPRSAx-»Q:(/i;) (compare Figs. |ll^d,e)). Since we consider lin- 
ear coupling, BPRSA5^_^Q,(A;) = — BPRSA^_^q(/c) as discussed in Section 3.2 
and illustrated in Figs. IH^d,e). There is no advantage over CCFx,y(fc) which 
looks very similar in this example. 

4-2 Nonlinear relation 

The response of the BPRSA to nonlinearly coupled trigger and target signals 
strongly depends on the type of the coupling. The most simple nonlinear 
coupling is the absolute value. Let us assume a sinusoidal trigger signal X 
without noise and a target signal Y that only contains the absolute value of 
X, yielding a frequency doubling. When calculating the BPRSA all oscillations 

cancel out and BPRSAx .y(A;) = BPRSAy ,x{k) = 0. In the presence 

of additional l//-noise the BPRSA will basically show features of the noise 
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Fig. 4. Samples of the noise series X (a, pure noise), Y (b, generated from X by 
Eq. ([II]) with ci = 0.2), and Z (c, ci = -0.1). The HF band (/ G [0.25,0.35] 
reciprocal sampling units) is used for the bandpass filtering, and the total length of 
the data is iV = 16384. BPRSA results for a = y (d) and a = Z (e): BPRSA^^^ 
(black sohd lines), BPRSA£^^ (red solid hues), - BPRSA^^^ ( green triangles), 
— BPRSA^^ (blue circles) are shown. The points are connected for visual reasons 
only; all values are dimensionless. 

and possibly finite size effects. The same holds for similar nonlinear coupling, 
e.g., raising to an even power. On the other hand, this elimination of higher 
harmonics might be an advantage if one wants to clarify a complex relationship 
between two unknown signals. 

Now, we study nonlinear coupling without frequency doubling. Three simple 
oscillating series are defined by 

= sin(27r/z), y^ = {xif, Zi = sgn{xi)\xi\^''^ . (12) 

and illustrated in Figs. [5](a-c). The large power of 9 has been chosen for visual 
reasons only; it enhances the differences as does the absence of noise. The cross- 
correlation analysis (see Figs. [5]^f,g)) cannot distinguish (i) the cases X —^Y 
and X — > Z as well as (ii) both possible analysis directions. Studying only the 
cross-correlation function could thus lead to the false conclusion of equivalently 
related signals Y and Z. BPRSA, on the other hand, can clearly distinguish 
the four cases except for BPRSAy^x(^) = BPRSA^^x(^)- However, one 
has to keep in mind that the shape of the BPRSA curve needs not be the 
same as the original target signal (compare Figs. [5]^b,d)). A presence of noise 
might disturb the BPRSA signal, making the identification of characteristics 
in trigger and target signal more difficult, depending on the signal to noise 
ratio. 
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Fig. 5. (a) Sinusoidal signal X and nonlinearly coupled signals (b) Y and (c) Z 
according to Eq. 1^. (d) BPRSA5^_^y (black solid lines), BPRSA^_^^ (red solid 
lines), -BPRSA^^y (green triangles), - BPRSA^__^^ (blue circles), (e) BPRSA 
for Z instead of Y accordingly; panel (f) shows CCFx,y = CCFy^x and (g) 
CCFx,z = CCFz,x (black on red solid lines) 

4-3 Influence of Trends in the signal 

Now let X and Y be two independent l//-noise signals with zero mean and 
unit variance generated by Fourier filtering. Furthermore, we add to both 
signals a periodic component Asin{27rfi). Moreover, non-stationarities are 
introduced by adding piecewise linear trends as follows. We start with some 
initial value for the slope ai and the initial offset Oq. At random positions, 
the offset and the slope are changed randomly within a previously defined 
range; the trends added to X and Y are independent (see Figs. [6](a,b)). For 
comparison we define a third signal Z that equals Y without trends (Fig. 
He)). 

Trends in the trigger signal will hardly affect the identification of the anchor 
points, because the anchor criteria defined in Eq. ([1]) is only based on local 
fluctuations. Note, that this might be different when using a more sophisti- 
cated boolean anchor function as discussed earlier (compare BPRSA directions 
r ^ X and Z ^ X in Fig. [S](d,e)). 
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Fig. 6. Sinusoids with frequency / = 1/30, amplitude A = 0.5, and normalized ad- 
ditional 1/f noise with /3 = 1.0; length N = 16384. In (a,b), different partial trends 
of random offset, slope and duration were added; (c) is the same as (b) but with- 
out trends, (d) BPRSA results for Y: BPRSA5^^y (black sohd hues), BPRSA^^^ 
(red solid lines), — BPRSA^^y (green triangles), — BPRSAy^j^ (blue circles), (e) 
BPRSA for Z replacing Y accordingly. Panel (f) shows CCFx,y (black), CCFy^x 
(red) and (g) CCFx,z (black), CCFz,x (red) accordingly. The points are connected 
for visual reasons only; all values are dimensionless. 

On the other hand, the influence of trends in the target signal cannot be 
neglected (see Fig. [6](e)). In case of a significant global trend in the tar- 
get signal, e.g., more decreasing parts than increasing parts, the global 
trend will be present in the BPRSA curve, although it is diminished. Note, 

that due to trends which do not cancel out completely BPRSAj^ >y(^) 7^ 

— BPRSA^ >y(^) in general (compare solid lines and triangles in Figs. 

[6]^d,e)). When the BPRSA shows no trend at all, the target signal is either 
characterized by no trends or the duration and slopes of increasing and de- 
creasing trends cancel out. 

As an implication of the different influences of trends in the trigger and tar- 
get signal one can identify which signal is disturbed by trends by comparing 
the BPRSA for opposite trigger-target directions {X ^ Y, Y ^ X). This 
is inherently impossible with cross-correlation analysis since the algorithm 
does not distinguish between both signals. Besides, trends are harmful for the 
definition of a global mean and thus disturb the standard cross-correlation 
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analysis. Therefore, its results may suggest a wrong correlation behavior. In 
Figs. [6](f),(g) one finds, by chance, anti-correlated behavior although the sig- 
nals themselves, i.e., the sinusoids, are strongly positively correlated. For the 
same reason a normalized BPRSA as defined in Eq. ([H]) cannot be applied 
here. Of course, in this simple example the use of the local cross-correlation 
function, which is based on local means rather than on a global mean, might 
help to remove the infiuence of the trends. 



5 Summary and Outlook 

In summary we have shown that the BPRSA method has several advan- 
tages compared with conventional cross-correlation analysis in the detection 
of quasi-periodicities in noisy non-stationary data with oscillations of finite 
coherence time. The method allows the analysis of the inter-relationship be- 
tween two signals, in particular effects in one signal triggered by events in 
another signal. 

This capability can be useful for the study of the inter-relation between res- 
piration, heart rate and blood pressure, i.e., the cardiovascular regulation, 
which is an important topic in human physiology. Cardiovascular functions are 
controlled by the tone of the sympathetic and parasympathetic (autonomic) 
nervous system that is infiuenced by the barorefiex, a homeostatic regulation 
that maintains a 'stable' blood pressure. An elevated blood pressure refiexively 
causes the blood pressure to decrease and vice versa. It is controlled through 
several stretch sensitive mechanoreceptors (baroreceptorslHI. It is believed that 
cardiovascular illnesses disturb the barorefiex. Related parameters might thus 
improve currently used predictors. Hence, the detection of quasi-periodicities 
refiecting regulation processes of the autonomic cardiac nervous system coin- 
ciding with increases or decreases of blood pressure in long records of human 
heart rate is of high clinical relevance. Autonomic dysfunction is closely related 
to cardiac mortality and susceptibility to life-threatening arrhythmic events 
[IB] . The assessment of heart rate variability by the PRSA based deceleration 
capacity (DC) parameter pUfTT] was shown to be superior to spectral param- 
eters proposed earlier for risk prediction pJ3j. BPRSA seems to be promising 
for the definition of an advanced risk predictor that respects the coupling of 



Activation of the baroreceptor results in an inhibition of sympathetic compo- 
nents and activation of parasympathetic or vagal components. Due to an initially 
elevated blood pressure activated baroreceptors tend to decrease cardiac output via 
a decrease in contractility resulting in a lower heart rate and finally in a decrease in 
blood pressure. A low blood pressure level relaxes the mechanoreceptor and stops 
the sympathetic inhibition and results in an increased contractility, heart rate and 
blood pressure. 
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heart rate variability and blood pressure. 

Further possible applications of BPRSA in biology and physiology include 
rhythmic motions of limbs in walking, muscle contractions, rhythms underly- 
ing the release of hormones that regulate growth and metabolism, periodicities 
in gene expression, membrane potential oscillations, oscillations in neuronal 
signals, and circadian rhythms ^lii2j . We believe that the range of suitable 
applications for the BPRSA method also includes quasi-periodic geophysical 
data, e.g., the El-Nino phenomenon, sunspot numbers, and ice age periods 
[3]. In addition, the analysis of complex elastic wave patterns to study seis- 
mic events or to determine material properties of granular matter might be 
improved by BPRSA. The study of non-stationary quasi-periodic complex 
waveforms is also a common task in the analysis and recognition of speech or 
music. 
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